#library(ncdf4)
library(ggplot2)
require(data.table)
library(dplyr)
library(plyr)
library(forecast)

#change to full local path to the replication data
wd1 <- "C:\\Users\\YFan\\OneDrive - NORCE\\NFood_replication\\"

#load intermediate Rdata prepared by "Data.Fig.S9.R"
load(paste0(wd1,"/data/RData/Fig.S9.GAM-effect-wtm.Rdata"))
gam.predict.ma <- gam.effect.wm[, lapply(.SD, ma, order=5), by =.(Scenario,Crop,Irrig,Level)] #by =.(Scenario,Crop,Irrig)] #
gam.predict.ma <- gam.predict.ma[, lapply(.SD, as.numeric), by =.(Scenario,Crop,Irrig,Level)] #by =.(Scenario,Crop,Irrig)] #
gam.predict.ma[,Year:=as.integer(Year)]
gam.predict.ma <- gam.predict.ma[!is.na(Year)]


#load dif.8crop.wm prepared by Data.Fig3-4.R
load(paste0(wd1,"/data/RData/Fig3-4.globalmean-yieldchanges.Rdata")) #including all grids with crop failures to calc modeled total change in yield under different scenarios (area>1000ha)
dif.8crop.wm.ma <- dif.8crop.wm[, lapply(.SD, ma, order=5), by =.(Scenario,Crop,Irrig)] #by =.(Scenario)] #
dif.8crop.wm.ma <- dif.8crop.wm.ma[, lapply(.SD, as.numeric), by =.(Scenario,Crop,Irrig)] #by =.(Scenario)] #
dif.8crop.wm.ma[,Year:=as.integer(Year)]
dif.8crop.wm.ma <- dif.8crop.wm.ma[!is.na(Year)]

eff.ma.pct <- gam.predict.ma[dif.8crop.wm.ma,
                                    .( Level, T.eff=(exp(T.log)-1)*100, #R.eff=(exp(R.log)-1)*100, M.eff=(exp(M.log)-1)*100,
                                       RD.eff=(exp(RD.log)-1)*100, RI.eff=(exp(RI.log)-1)*100, P.eff=(exp(P.log)-1)*100, RH.eff=(exp(RH.log)-1)*100,
                                       Tot.eff=(exp(Tot.log)-1)*100, co2.eff=(exp(co2.log)-1)*100, #luc.eff=(exp(luc.log)-1)*100, 
                                       Y.mod.dif=(exp(i.Y.mod.dif)-1)*100, #Y.mod.dif=(yield-yield0)/yield0*100, #differene of wtm yield is not the same as wtm of per-grid difference
                                       Tot.eff2=(exp(Tot.log2)-1)*100, Y.mod.dif2=(exp(i.Y.mod.dif2)-1)*100, #climate+CO2 for ER
                                       yield=i.yield, yield1=i.yield1, yield0=i.yield0, area=i.area 
                                    ), by=.EACHI, #on=.(Scenario,Year)]
                                    on=.(Scenario,Crop,Irrig,Year)]
eff.ma.pct <- eff.ma.pct[!is.na(Level)]
#aggregate per crop effect to global considering base yield difference between crop types
eff.tot.ma <- rbind(eff.ma.pct[Scenario=="RCP45 - RCP85",
                                     .(T.eff=weighted.mean(T.eff, yield1*area, na.rm = T), #R.eff=weighted.mean(R.eff, yield1*area, na.rm = T), M.eff=weighted.mean(M.eff, yield1*area, na.rm = T),
                                       RD.eff=weighted.mean(RD.eff, yield1*area, na.rm = T),RI.eff=weighted.mean(RI.eff, yield1*area, na.rm = T),
                                       P.eff=weighted.mean(P.eff, yield1*area, na.rm = T), RH.eff=weighted.mean(RH.eff, yield1*area, na.rm = T),
                                       Tot.eff=weighted.mean(Tot.eff, yield1*area, na.rm = T), Y.mod.dif=weighted.mean(Y.mod.dif, yield1*area, na.rm = T), 
                                       co2.eff=weighted.mean(co2.eff, yield0*area, na.rm = T), #luc.eff=weighted.mean(luc.eff, yield*area, na.rm = T), #LUC has differnet base yield from RCP45_85LU
                                       Tot.eff2=weighted.mean(Tot.eff2, yield0*area, na.rm = T), Y.mod.dif2=weighted.mean(Y.mod.dif2, yield0*area, na.rm = T),
                                       #Tot.eff3=weighted.mean(Tot.eff3, yield0*area, na.rm = T), Y.mod.dif3=weighted.mean(Y.mod.dif3, yield0*area, na.rm = T), 
                                       yield1=weighted.mean(yield1, area), yield=weighted.mean(yield, area), #yield2=weighted.mean(yield2, area2), area2=sum(area2, na.rm=T),
                                       yield0=weighted.mean(yield0, area), area=sum(area, na.rm=T)),
                                     by =.(Scenario,Level,Year)],
                    eff.ma.pct[Scenario!="RCP45 - RCP85",
                                     .(T.eff=weighted.mean(T.eff, yield0*area, na.rm = T), #R.eff=weighted.mean(R.eff, yield0*area, na.rm = T),M.eff=weighted.mean(M.eff, yield0*area, na.rm = T),
                                       RD.eff=weighted.mean(RD.eff, yield0*area, na.rm = T),RI.eff=weighted.mean(RI.eff, yield0*area, na.rm = T),
                                       P.eff=weighted.mean(P.eff, yield0*area, na.rm = T),RH.eff=weighted.mean(RH.eff, yield0*area, na.rm = T),
                                       Tot.eff=weighted.mean(Tot.eff, yield0*area, na.rm = T), Y.mod.dif=weighted.mean(Y.mod.dif, yield0*area, na.rm = T),
                                       yield0=weighted.mean(yield0, area), yield=weighted.mean(yield, area),
                                       area=sum(area, na.rm=T)),
                                     by =.(Scenario,Level,Year)], fill=TRUE)


eff.tot <- copy(eff.tot.ma)
eff.tot <- eff.tot[Scenario=="RCP45 - RCP85", Tot.eff:=Tot.eff2] #climate+CO2 only
eff.tot <- eff.tot[, Y.mod.dif:=(yield-yield0)/yield0*100] #relative change of global average yield (yield=RCP45_85LU,SAI,MSB,CCT, yield0=RCP85)

#compared tot.eff and Y.mod.dif
eff.tot[Year >=2075 & Year<=2099 & Level=="Mean", mean(Tot.eff), by=.(Scenario)]
eff.tot[Year >=2075 & Year<=2099 & Level=="Mean", mean(co2.eff), by=.(Scenario)]
eff.tot[Year >=2075 & Year<=2099 & Level=="Mean", mean(Y.mod.dif), by=.(Scenario)] 

#% effect on global yield
eff.tot.t <- eff.tot[, data.table(t(.SD), keep.rownames=TRUE), .SDcols=-"Level", by=.(Year,Scenario)]
#names(eff.tot.t)[3:6] <- c("Variable", "Mean", "Conf.L", "Conf.U")
names(eff.tot.t)[3:4] <- c("Variable", "Mean")



# library(RColorBrewer)
# cbPalette <- brewer.pal(n =12, "Paired") # ,"Dark2")
# #grays <- c("#D9D9D9", "#BDBDBD", "#969696", "#737373", "#525252", "#252525", "#000000")
# varcolor <- c(cbPalette[c(10,7,8,1,2,3)], "#000000", "#737373")
# barplot(height=rep(1,length(varcolor)),col=varcolor) #display colors

#====final colors for variables and scenarios
varcolor <- c("#E64B35FF","#FF7F00" ,"#FDBF6F", "#A6CEE3", "#1F78B4", "#00A087FF", "#000000", "#737373")
barplot(height=rep(1,length(varcolor)),col=varcolor) #display colors


linetype <- c(1,1,1, 1,1,1,1, 1)
size <- c(1,1,1, 1, 1, 1, 1, 1)*1
override.linetype <- 1
override.size <- 1.1
casename = c(expression('('*bold(a)*') SAI - RCP85'), expression('('*bold(b)*') MSB - RCP85'), expression('('*bold(c)*') CCT - RCP85'), 
             expression('('*bold(d)*') Emissions reduction')) #expression('('*bold(i)*') RCP45 - RCP85'))
cases = c("SAI - RCP85", "MSB - RCP85", "CCT - RCP85", "RCP45 - RCP85")
vars = c("T.eff", "RD.eff","RI.eff", "P.eff", "RH.eff", "co2.eff", "Tot.eff", "Y.mod.dif")
varname <- c(expression(atop(atop(textstyle('T'), textstyle('effect')), ' ')), #expression(atop("Diffuse \n radiation", 'effect')),
             expression(atop(atop(textstyle('RD'), textstyle('effect')),' ')),
             expression(atop(atop(textstyle('RI'),textstyle('effect')),' ')),
             expression(atop(atop(textstyle('P'), textstyle('effect')),  ' ')),
             expression(atop(atop(textstyle('RH'), textstyle('effect')),  ' ')),
             expression(atop(atop(textstyle('CO'[2]), textstyle('effect')), ' ')),
             expression(atop(atop(textstyle('Total'), textstyle('(MLR)')), ' ')),
             expression(atop(atop(textstyle('Total'), textstyle('(CLM5)')), ' ')))

p0 <-  eff.tot.t[Variable%in%vars] %>%
  mutate(Var = factor(Variable, levels=vars)) %>%
  mutate(Scenario = factor(Scenario, levels=cases, labels=casename)) %>%
  ggplot(aes(Year, Mean, group=Variable)) +
  geom_line(aes(color=Var, linetype=Var, size=Var)) + #use Case as group and set override values in manual scale
  geom_hline(aes(yintercept=0), color="gray") +
  #geom_linerange(aes(x=Year, ymin=Conf.L, ymax=Conf.U, color=Var), alpha=0.3, show.legend=FALSE) +
  #change linetype, color, size mannually
  scale_linetype_manual(values=linetype, guide = FALSE)+
  scale_size_manual(values=size, guide = FALSE) +
  scale_colour_manual(name="", labels=varname, values=varcolor)+ # cbPalette[c(1:4,4,5,5)]) + #cbPalette[c(1:5,5,6,6)]) +
  #use free scales for each panel
  facet_wrap(~Scenario, nrow=2, scales='fixed', labeller= label_parsed) +
  scale_x_continuous(name=NULL, limits = c(2020,2099), breaks =c(2025,2050,2075,2099)) +
  scale_y_continuous("Change in global yield (%)")+# limits = c(NA,40))+#, breaks = c(-10,0,10,20)) + #, sec.axis = dup_axis())
  guides(colour = guide_legend(keywidth =1,keyheight =0.8, nrow=1, byrow=T, label.position = "right", label.hjust = 0, label.vjust = 0,
                               override.aes = list(size = override.size, linetype = override.linetype ))) #merge the legends for color and size
p0 <- p0 + theme_bw() + #theme_minimal() + #
  ggtitle("") + theme(plot.title = element_text(size=14, face="bold", hjust = 0.5)) +
  theme(panel.spacing.x = unit(1, "lines")) +  #facet spacing
  theme(strip.text = element_text(size = 13)) +
  theme(legend.text = element_text(size=12), legend.title = element_blank(), legend.position = "bottom") +
  theme(axis.text = element_text(size=13), axis.title=element_text(size=13))
# set transparency for both plot and panel background
p0 + theme(strip.background = element_blank(),
           panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
           panel.background = element_rect(fill = "transparent",colour = NA),
           plot.background = element_rect(fill = "transparent",colour = NA))


ggsave(paste0(wd1,"/figures/supplementary/Fig.S9.svg"), device ="svg",  #vertical plots
       units = "in", width = 7, height = 6.5)




